Incremental response of granular materials: DEM results 
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Abstract. We systematically investigate the incremental response of various equilibrium states of dense 2D model granular 
materials, along the biaxial compression path (0\\ < 022, 012 = 0). Stress increments are applied in arbitrary directions in 3- 
dimensional stress space (oi i , O22, OI2). In states with stable contact networks we compute the stiffness matrix and the elastic 
moduli, and separate elastic and irreversible strains in the range in which the latter are homogeneous functions of degree 
one of stress increments. Without principal stress axis rotation, the response abides by elastoplasticity with a Mohr-Coulomb 
criterion and a non-associated flow rule. However a nonelastic shear strain is also observed for increments of G\2, and shear 
and in-plane responses couple. This behavior correlates to the distribution of friction mobilization and sliding at contacts. 
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INTRODUCTION 

Although the mechanical behavior of solidlike granular 
materials under quasistatic loading conditions is often 
modeled as elastoplastic at the continuum level [1, 2], 
there are still few studies addressing the microscopic 
origins of such a behavior by discrete, grain-level sim- 
ulation [3-6]. To assess the applicability of elastoplas- 
tic laws, one needs to investigate the response to small 
stress or strain increments, superimposed in various di- 
rections on an equilibrium state. One essential motiva- 
tion for such studies is the prediction of shear band for- 
mation, for which such criteria as the Rudnicki-Rice [7] 
condition involve the incremental response. In particu- 
lar, localization is crucially sensitive to the response to 
stress increments with rotation of principal axes, as when 
some simple shear is superimposed on a biaxial compres- 
sion [8]. The present study addresses this issue for the 
simplest model material, an assembly of disks in 2 di- 
mensions, for which the response to load increments in 
all 3 dimensions of stress space is computed at various 
points along a biaxial loading path. 



MODEL MATERIAL AND METHODS 

Our simulation samples comprise 5600 disks enclosed 
in a periodic rectangular cell. The diameter distribution 
is uniform between O.ld and 1.3d. We use a simple, 
frictional-elastic contact model, involving (constant) nor- 
mal contact stiffness Kn, tangential contact stiffness Kj 
(here we set Kj — Km) and a friction coefficient, fi, set 
to 0.3. The normal (elastic) contact force is Fn = K^h 
where /; is the interpenetration of contacting disks (which 
models surface deflection). The tangential force Fj re- 



lates to the elastic part 8 of the tangential relative dis- 
placement, as Ft — Kt8, and is incrementally computed 
to enforce the Coulomb condition \Ft\ < IJ.Fn- Some vis- 
cous damping is also introduced, which proves irrelevant 
to the material behavior for low enough strain rates. 

We focus here on dense samples, which are initially 
assembled without friction, under an isotropic pressure 
P. The initial state is thus characterized by an isotropic 
fabric and a large coordination number (close to 4). The 
dimensionless stiffness parameter k = K^/P sets the 
scale of contact deflections, as h/d °^ K . We choose 
value K — 10 4 in most simulations. 

Deformations of the simulation cell, i.e. macroscopic 
strains, are controlled, or vary in response to applied 
stresses. This is achieved with specific implementa- 
tions of Parrinello-Rahman and Lees-Edwards tech- 
niques (first developed for molecular systems [9]), as ex- 
plained in Ref. [10]. Stresses are given by the classical 
Love formula. In the biaxial compression test, the de- 
formable cell remains rectangular, its edges parallel to 
the principal stress directions. Principal stress value <j\ 
(the lateral stress) is kept equal to P, while 02 (the ax- 
ial stress) increases in response to strain £2, which grows 
at a controlled rate (compressive stresses and shrinking 
strains are positive). As indicated in Fig. 1, the com- 
pression test is stopped at different stages and the sam- 
ple is equilibrated at constant stresses. This entails slight 
creep strain increments, which remain quite small (of or- 
der 10~ 6 ), until equilibrium conditions are satisfied with 
good accuracy (the tolerance is 10~ 4 in units of P, dP, 
and d 2 P for stresses, forces and moments, respectively). 
In those well-equilibrated intermediate states, hereafter 
referred to as investigation points, we first compute elas- 
tic moduli. To do so, we use the stiffness matrix associ- 
ated to the contact network, as in Ref. [11]. It is conve- 
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FIGURE 1. Deviator stress versus axial strain curve, and 
location of 2 investigation points for incremental response 



nient to denote stresses and strains as 3-vectors (as 5(7, 
5e) with 5(73 = v2oi2> while notations 5<7i, 5<72 keep 
the same meaning (and similarly for 5e). Due to sym- 
metry about the principal axes there are four independent 
elastic moduli, which satisfy: 



5cr = C-5e £ , with C 
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superscripts E recalling that strains are purely elastic. 
The incremental response for various load directions 
is then computed, for different stress increments. We 
choose 5(7 values on a sphere in 3-space, centered at the 
origin, of radius 2\/2 x 10~ 3 P. Such increments are ap- 
plied, and then multiplied by integer factors 2, 3... up 
to 12, in order to record the influence of both their di- 
rection and their amplitude. The calculations are fully 
stress-controlled, with variations of all 3 strain compo- 
nents. Once a new, pertubed equilibrium is reached, 5e 
is measured, from which the elastic part de E = C ■ 5(7 
is subtracted, defining the irreversible strain increment, 
which we denote with superscript P (for "plastic"). In- 
vestigation points with Oil <5\ = 1.2, 1.4, 1.6 and 1.8 
were studied 

Before presenting the results in the next section, let 
us recall that as a consequence of the assembling pro- 
cess, the investigated states possess a large coordination 
number, and (see Fig. 1) are within the range of strain 
("regime I"), along the biaxial loading curve, that is dom- 
inated by contact deformation [11, 12]. This means that 
the contact network does not break apart, and that the 
irreversible strains are due to sliding at contacts where 
the Coulomb limit is reached. As a consequence, macro- 
scopic strains, on changing confining stresses or stiff- 
ness constants, scale as fc [12]. For larger deviators, 
or in poorly coordinated samples (which might be very 
dense nevertheless [11]), the macroscopic strains stem 
from network ruptures and rearrangements ("regime II") 
and the incremental behavior might differ. 
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FIGURE 2. Stress increments in principal axis plane 



INCREMENTAL RESPONSE 
No rotation of principal axes 

We first investigate the response to stress increments 
lying in the plane of principal stresses (i.e., 5(73 = 0). 
For each investigation point along the curve of Fig. 1, 12 
different orientations of 5(7 in this plane are tested, as 
shown in Fig. 2, with 12 different amplitudes (as speci- 
fied before). In order to assess the relevance of classical 
plasticity models for the material studied here we focus 
on the following three aspects: (i) the existence of a flow 
rule dictating the direction of 5e p ; (ii) at equal ampli- 
tude |5ff|, the linear dependence of amplitude |5e p | on 
the positive part [5(7]+ of 5(7 = Nc • 5(7, where Nc is 
the outer normal to some yield criterion in stress space; 
(iii) the same linear dependence for varying stress incre- 
ment amplitudes. The existence of a plastic flow rule is 
a sharp feature arising from incremental tests, as shown 
in Fig. 3, corresponding to an investigation point with 
O2/01 = Elastic strain increments 5ef and 5&f are 
disposed along as many directions as the stress incre- 
ments in Fig. 2, while plastic strain increments (5ef and 
5&f) clearly align along a unique direction, consistently 
with the flow rule. The same features are observed for all 
investigation points. 

We discuss point (ii) of our list by referring to Fig. 4 
in which |5e p | is plotted versus the angle a between 
principal axis 1 and increment 5 a, at constant amplitude 
1 5(7 1 . In the framework of classical plasticity these values 
should fit to the positive part of a cosine function reach- 
ing its maximum in tnormal to the yield criterion. Fit- 
ting theoretical curves to data allows to estimate the an- 
gle CCnyc characterising the normal Nc to the yield crite- 
rion [4] and the maximal amplitude de^ AX of the plastic 
strain increment. Notably, for all investigation points, the 
normal Nc, oriented at angle a^yc, is consistently very 
nearly orthogonal to the current stress direction (7i,(72 
(oriented at angle o,ld in stress space). This suggests 
that the yield criterion might be defined by the Coulomb 
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FIGURE 3. Elastic and anelastic parts of response to stress 
increments marked (Oa, Ob, 01) in Fig. 2 




FIGURE 4. Amplitude \8e F 



orientation a of 8 a for 



constant amplitude \So\ = 3.394- 10~ 2 (O2M = 1-4). 



condition of a constant ratio a%jo\. Since OCpfd 7^ (%nyc 
the plastic flow direction differs from the normal No as 
in nonassociated elastoplasticity . As to point (iii), it is 



0.8 
0.7 



0.61 

a, 

TW 0.5- 

0.3 
0.2 
0.1 
0.0 



+ CT2/CT1 = 1-2 
x CT2/C1 = 1.4 
© CT2/C1 = 1-6 
■»• O2/C1 = 1 



0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035 

FIGURE 5. Nonelastic strain amplitude vs. defined 
with normal to criterion identified in Fig. 4. 

checked in Fig. 5, from which the following plastic mod- 
uli Cp (in units of Kpf) are measured: Cp =??, ??, ??, ?? 
corresponding, respectively, to G2/01 =1.2, 1.4, 1.6 and 
1.8. 



General case 

If elastoplasticity applies - which seems to be the 
case for 5(7 in the plane of the principal stress direc- 
tions - then a small load increment in the third direction, 
5(73 7^ 0, 8(J\ = 5(72 = should entail a purely elastic 
response. Fig. 6 contradicts this prediction, as a nonelas- 
tic shear strain 8e^ immediately appears, which increase 
proportionnally to shear stress | cTi 2 1 . Coefficients can be 
slightly different for positive and negative 5(7i2 because 
of finite sample size effects. Like in-plane increments, 
such 5(7, if extremely small, yield a nonelastic response 
that is slightly sublinear in their amplitude, but a plas- 
tic modulus can be identified for 5(73 /P of order 10~ 2 . 
Out-of-plane increments 5(7 also entail plastic strains 
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FIGURE 6. Total, elastic, nonelastic shear strains as func- 
tions of applied shear stress to state with 02/ d\ = 1.8. Plastic 
modulus is close to 3Kpj (resp. for Son > (5(7i2 < 0). 

5ef, deJ? , which are still related by the same flow rule 
as previously identified for in-plane loads (5(73 = 0). 
Fig. 7 gathers results both from 16 load directions for 

which 5(73 = ±-y/ 5(7 X 2 + 5(72' as wei ^ as simple shear 
increments (5(7i = 5(72 =0) with both signs of 5(73. 
Quite surprisingly, the latter also produce a nonelastic re- 
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FIGURE 7. Analog of Fig. 3 in state with (7 2 /(7i = 1-8, for 
out-of-plane So. Big red dots correspond to 8<T\ = 802 = 0. 
Elastic strains (bottom right) are comparatively smaller. 



ponse in the plane of principal stresses. We thus observe 



that both the irreversible strains and the stress incre- 
ments causing them span two-dimensional spaces, with 
one in-plane and one out-of-plane direction, and that the 
response couples both directions. To be complete, we 
should then specify how 8e F depends on 8d for all load 
increments. Although we are still investigating this issue, 
some preliminary attempts at superposition of responses 
to shear and to in-plane stress increments are encourag- 
ing, as shown by Fig. 8. Upon superimposing the pre- 
cis A 
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FIGURE 8. Predicted, with procedure defined in text, versus 
observed 8e2 for combined loads (o^/ci = 1.8). 

viously identified responses to (in-plane) 86 = Nc • 8d 
and to 1 5(73 1 in simple shear, Fig. 8 shows that the pre- 
dicted values are fairly close to the measured ones. 

MICROSCOPIC ASPECTS 

The macroscopic nonelastic is due to plastic sliding in 
some contacts. While the distribution of contact orien- 
tations (fabric) is still moderately anisotropic in the in- 
vestigated states, the sliding contact fabric (Fig. 9) has a 
much stronger angular dependence. Such a distribution 




FIGURE 9. Left: sliding contact orientational distribution 
(major principal axis vertical on the plot), normalized such that 
its angular average is a coordination number. Diameter of circle 
is 1.1, global coordination is 3.5. Right: angular distribution 
of amplitude of sliding relative displacement in contacts in 
response to in-plane load increment, normalized by plastic 
strain. &ija\ = 1.8. 



is observed for 85 > 0, while the population of slid- 
ing contacts virtually vanishes on applying 86 < and 
5(73 =0. The sliding contact fabric depends on both 5(73 
and 5(7 in general. A nonzero 5(73 breaks its symmetry. 
The angular distribution of sliding displacements at con- 
tacts (Fig. 9), albeit different, is also strongly anisotropic 
and shows similar sensitivity to the direction of 5(7. Fi- 
nally, stress increments for which 5(7 is proportional to a 
(the neutral direction), entail no sliding, as contact forces 
tend to increase proportionnally to their previous value. 

PERSPECTIVES 

The essential finding of the present study, which still re- 
mains to be systematized and calls for more thorough mi- 
cromechanical investigations, is the correspondence be- 
tween 2D stress increments orthogonal to the currrent 
stress level and nonelastic strains belonging to a 2D 
space. In the near future we plan to formulate it as a com- 
plete constitutive incremental law, to relate it to micro- 
scopic phenomena and to use it in localization criteria. 
The incremental response in systems with gradually re- 
arranging contact networks ("regime II", associated with 
microscopic instabilities) should also be investigated. 
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